home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dhseqr.f < prev    next >
Text File  |  1996-07-19  |  16KB  |  455 lines
  1.       SUBROUTINE DHSEQR( JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z,
  2.      $                   LDZ, WORK, LWORK, INFO )
  3. *
  4. *  -- LAPACK routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       CHARACTER          COMPZ, JOB
  11.       INTEGER            IHI, ILO, INFO, LDH, LDZ, LWORK, N
  12. *     ..
  13. *     .. Array Arguments ..
  14.       DOUBLE PRECISION   H( LDH, * ), WI( * ), WORK( * ), WR( * ),
  15.      $                   Z( LDZ, * )
  16. *     ..
  17. *
  18. *  Purpose
  19. *  =======
  20. *
  21. *  DHSEQR computes the eigenvalues of a real upper Hessenberg matrix H
  22. *  and, optionally, the matrices T and Z from the Schur decomposition
  23. *  H = Z T Z**T, where T is an upper quasi-triangular matrix (the Schur
  24. *  form), and Z is the orthogonal matrix of Schur vectors.
  25. *
  26. *  Optionally Z may be postmultiplied into an input orthogonal matrix Q,
  27. *  so that this routine can give the Schur factorization of a matrix A
  28. *  which has been reduced to the Hessenberg form H by the orthogonal
  29. *  matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.
  30. *
  31. *  Arguments
  32. *  =========
  33. *
  34. *  JOB     (input) CHARACTER*1
  35. *          = 'E':  compute eigenvalues only;
  36. *          = 'S':  compute eigenvalues and the Schur form T.
  37. *
  38. *  COMPZ   (input) CHARACTER*1
  39. *          = 'N':  no Schur vectors are computed;
  40. *          = 'I':  Z is initialized to the unit matrix and the matrix Z
  41. *                  of Schur vectors of H is returned;
  42. *          = 'V':  Z must contain an orthogonal matrix Q on entry, and
  43. *                  the product Q*Z is returned.
  44. *
  45. *  N       (input) INTEGER
  46. *          The order of the matrix H.  N >= 0.
  47. *
  48. *  ILO     (input) INTEGER
  49. *  IHI     (input) INTEGER
  50. *          It is assumed that H is already upper triangular in rows
  51. *          and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
  52. *          set by a previous call to DGEBAL, and then passed to SGEHRD
  53. *          when the matrix output by DGEBAL is reduced to Hessenberg
  54. *          form. Otherwise ILO and IHI should be set to 1 and N
  55. *          respectively.
  56. *          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
  57. *
  58. *  H       (input/output) DOUBLE PRECISION array, dimension (LDH,N)
  59. *          On entry, the upper Hessenberg matrix H.
  60. *          On exit, if JOB = 'S', H contains the upper quasi-triangular
  61. *          matrix T from the Schur decomposition (the Schur form);
  62. *          2-by-2 diagonal blocks (corresponding to complex conjugate
  63. *          pairs of eigenvalues) are returned in standard form, with
  64. *          H(i,i) = H(i+1,i+1) and H(i+1,i)*H(i,i+1) < 0. If JOB = 'E',
  65. *          the contents of H are unspecified on exit.
  66. *
  67. *  LDH     (input) INTEGER
  68. *          The leading dimension of the array H. LDH >= max(1,N).
  69. *
  70. *  WR      (output) DOUBLE PRECISION array, dimension (N)
  71. *  WI      (output) DOUBLE PRECISION array, dimension (N)
  72. *          The real and imaginary parts, respectively, of the computed
  73. *          eigenvalues. If two eigenvalues are computed as a complex
  74. *          conjugate pair, they are stored in consecutive elements of
  75. *          WR and WI, say the i-th and (i+1)th, with WI(i) > 0 and
  76. *          WI(i+1) < 0. If JOB = 'S', the eigenvalues are stored in the
  77. *          same order as on the diagonal of the Schur form returned in
  78. *          H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2
  79. *          diagonal block, WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and
  80. *          WI(i+1) = -WI(i).
  81. *
  82. *  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
  83. *          If COMPZ = 'N': Z is not referenced.
  84. *          If COMPZ = 'I': on entry, Z need not be set, and on exit, Z
  85. *          contains the orthogonal matrix Z of the Schur vectors of H.
  86. *          If COMPZ = 'V': on entry Z must contain an N-by-N matrix Q,
  87. *          which is assumed to be equal to the unit matrix except for
  88. *          the submatrix Z(ILO:IHI,ILO:IHI); on exit Z contains Q*Z.
  89. *          Normally Q is the orthogonal matrix generated by DORGHR after
  90. *          the call to DGEHRD which formed the Hessenberg matrix H.
  91. *
  92. *  LDZ     (input) INTEGER
  93. *          The leading dimension of the array Z.
  94. *          LDZ >= max(1,N) if COMPZ = 'I' or 'V'; LDZ >= 1 otherwise.
  95. *
  96. *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
  97. *
  98. *  LWORK   (input) INTEGER
  99. *          This argument is currently redundant.
  100. *
  101. *  INFO    (output) INTEGER
  102. *          = 0:  successful exit
  103. *          < 0:  if INFO = -i, the i-th argument had an illegal value
  104. *          > 0:  if INFO = i, DHSEQR failed to compute all of the
  105. *                eigenvalues in a total of 30*(IHI-ILO+1) iterations;
  106. *                elements 1:ilo-1 and i+1:n of WR and WI contain those
  107. *                eigenvalues which have been successfully computed.
  108. *
  109. *  =====================================================================
  110. *
  111. *     .. Parameters ..
  112.       DOUBLE PRECISION   ZERO, ONE, TWO
  113.       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
  114.       DOUBLE PRECISION   CONST
  115.       PARAMETER          ( CONST = 1.5D+0 )
  116.       INTEGER            NSMAX, LDS
  117.       PARAMETER          ( NSMAX = 15, LDS = NSMAX )
  118. *     ..
  119. *     .. Local Scalars ..
  120.       LOGICAL            INITZ, WANTT, WANTZ
  121.       INTEGER            I, I1, I2, IERR, II, ITEMP, ITN, ITS, J, K, L,
  122.      $                   MAXB, NH, NR, NS, NV
  123.       DOUBLE PRECISION   ABSW, OVFL, SMLNUM, TAU, TEMP, TST1, ULP, UNFL
  124. *     ..
  125. *     .. Local Arrays ..
  126.       DOUBLE PRECISION   S( LDS, NSMAX ), V( NSMAX+1 ), VV( NSMAX+1 )
  127. *     ..
  128. *     .. External Functions ..
  129.       LOGICAL            LSAME
  130.       INTEGER            IDAMAX, ILAENV
  131.       DOUBLE PRECISION   DLAMCH, DLANHS, DLAPY2
  132.       EXTERNAL           LSAME, IDAMAX, ILAENV, DLAMCH, DLANHS, DLAPY2
  133. *     ..
  134. *     .. External Subroutines ..
  135.       EXTERNAL           DCOPY, DGEMV, DLABAD, DLACPY, DLAHQR, DLARFG,
  136.      $                   DLARFX, DLASET, DSCAL, XERBLA
  137. *     ..
  138. *     .. Intrinsic Functions ..
  139.       INTRINSIC          ABS, MAX, MIN
  140. *     ..
  141. *     .. Executable Statements ..
  142. *
  143. *     Decode and test the input parameters
  144. *
  145.       WANTT = LSAME( JOB, 'S' )
  146.       INITZ = LSAME( COMPZ, 'I' )
  147.       WANTZ = INITZ .OR. LSAME( COMPZ, 'V' )
  148. *
  149.       INFO = 0
  150.       IF( .NOT.LSAME( JOB, 'E' ) .AND. .NOT.WANTT ) THEN
  151.          INFO = -1
  152.       ELSE IF( .NOT.LSAME( COMPZ, 'N' ) .AND. .NOT.WANTZ ) THEN
  153.          INFO = -2
  154.       ELSE IF( N.LT.0 ) THEN
  155.          INFO = -3
  156.       ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
  157.          INFO = -4
  158.       ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
  159.          INFO = -5
  160.       ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
  161.          INFO = -7
  162.       ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. LDZ.LT.MAX( 1, N ) ) THEN
  163.          INFO = -11
  164.       END IF
  165.       IF( INFO.NE.0 ) THEN
  166.          CALL XERBLA( 'DHSEQR', -INFO )
  167.          RETURN
  168.       END IF
  169. *
  170. *     Initialize Z, if necessary
  171. *
  172.       IF( INITZ )
  173.      $   CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
  174. *
  175. *     Store the eigenvalues isolated by DGEBAL.
  176. *
  177.       DO 10 I = 1, ILO - 1
  178.          WR( I ) = H( I, I )
  179.          WI( I ) = ZERO
  180.    10 CONTINUE
  181.       DO 20 I = IHI + 1, N
  182.          WR( I ) = H( I, I )
  183.          WI( I ) = ZERO
  184.    20 CONTINUE
  185. *
  186. *     Quick return if possible.
  187. *
  188.       IF( N.EQ.0 )
  189.      $   RETURN
  190.       IF( ILO.EQ.IHI ) THEN
  191.          WR( ILO ) = H( ILO, ILO )
  192.          WI( ILO ) = ZERO
  193.          RETURN
  194.       END IF
  195. *
  196. *     Set rows and columns ILO to IHI to zero below the first
  197. *     subdiagonal.
  198. *
  199.       DO 40 J = ILO, IHI - 2
  200.          DO 30 I = J + 2, N
  201.             H( I, J ) = ZERO
  202.    30    CONTINUE
  203.    40 CONTINUE
  204.       NH = IHI - ILO + 1
  205. *
  206. *     Determine the order of the multi-shift QR algorithm to be used.
  207. *
  208.       NS = ILAENV( 4, 'DHSEQR', JOB // COMPZ, N, ILO, IHI, -1 )
  209.       MAXB = ILAENV( 8, 'DHSEQR', JOB // COMPZ, N, ILO, IHI, -1 )
  210.       IF( NS.LE.2 .OR. NS.GT.NH .OR. MAXB.GE.NH ) THEN
  211. *
  212. *        Use the standard double-shift algorithm
  213. *
  214.          CALL DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
  215.      $                IHI, Z, LDZ, INFO )
  216.          RETURN
  217.       END IF
  218.       MAXB = MAX( 3, MAXB )
  219.       NS = MIN( NS, MAXB, NSMAX )
  220. *
  221. *     Now 2 < NS <= MAXB < NH.
  222. *
  223. *     Set machine-dependent constants for the stopping criterion.
  224. *     If norm(H) <= sqrt(OVFL), overflow should not occur.
  225. *
  226.       UNFL = DLAMCH( 'Safe minimum' )
  227.       OVFL = ONE / UNFL
  228.       CALL DLABAD( UNFL, OVFL )
  229.       ULP = DLAMCH( 'Precision' )
  230.       SMLNUM = UNFL*( NH / ULP )
  231. *
  232. *     I1 and I2 are the indices of the first row and last column of H
  233. *     to which transformations must be applied. If eigenvalues only are
  234. *     being computed, I1 and I2 are set inside the main loop.
  235. *
  236.       IF( WANTT ) THEN
  237.          I1 = 1
  238.          I2 = N
  239.       END IF
  240. *
  241. *     ITN is the total number of multiple-shift QR iterations allowed.
  242. *
  243.       ITN = 30*NH
  244. *
  245. *     The main loop begins here. I is the loop index and decreases from
  246. *     IHI to ILO in steps of at most MAXB. Each iteration of the loop
  247. *     works with the active submatrix in rows and columns L to I.
  248. *     Eigenvalues I+1 to IHI have already converged. Either L = ILO or
  249. *     H(L,L-1) is negligible so that the matrix splits.
  250. *
  251.       I = IHI
  252.    50 CONTINUE
  253.       L = ILO
  254.       IF( I.LT.ILO )
  255.      $   GO TO 170
  256. *
  257. *     Perform multiple-shift QR iterations on rows and columns ILO to I
  258. *     until a submatrix of order at most MAXB splits off at the bottom
  259. *     because a subdiagonal element has become negligible.
  260. *
  261.       DO 150 ITS = 0, ITN
  262. *
  263. *        Look for a single small subdiagonal element.
  264. *
  265.          DO 60 K = I, L + 1, -1
  266.             TST1 = ABS( H( K-1, K-1 ) ) + ABS( H( K, K ) )
  267.             IF( TST1.EQ.ZERO )
  268.      $         TST1 = DLANHS( '1', I-L+1, H( L, L ), LDH, WORK )
  269.             IF( ABS( H( K, K-1 ) ).LE.MAX( ULP*TST1, SMLNUM ) )
  270.      $         GO TO 70
  271.    60    CONTINUE
  272.    70    CONTINUE
  273.          L = K
  274.          IF( L.GT.ILO ) THEN
  275. *
  276. *           H(L,L-1) is negligible.
  277. *
  278.             H( L, L-1 ) = ZERO
  279.          END IF
  280. *
  281. *        Exit from loop if a submatrix of order <= MAXB has split off.
  282. *
  283.          IF( L.GE.I-MAXB+1 )
  284.      $      GO TO 160
  285. *
  286. *        Now the active submatrix is in rows and columns L to I. If
  287. *        eigenvalues only are being computed, only the active submatrix
  288. *        need be transformed.
  289. *
  290.          IF( .NOT.WANTT ) THEN
  291.             I1 = L
  292.             I2 = I
  293.          END IF
  294. *
  295.          IF( ITS.EQ.20 .OR. ITS.EQ.30 ) THEN
  296. *
  297. *           Exceptional shifts.
  298. *
  299.             DO 80 II = I - NS + 1, I
  300.                WR( II ) = CONST*( ABS( H( II, II-1 ) )+
  301.      $                    ABS( H( II, II ) ) )
  302.                WI( II ) = ZERO
  303.    80       CONTINUE
  304.          ELSE
  305. *
  306. *           Use eigenvalues of trailing submatrix of order NS as shifts.
  307. *
  308.             CALL DLACPY( 'Full', NS, NS, H( I-NS+1, I-NS+1 ), LDH, S,
  309.      $                   LDS )
  310.             CALL DLAHQR( .FALSE., .FALSE., NS, 1, NS, S, LDS,
  311.      $                   WR( I-NS+1 ), WI( I-NS+1 ), 1, NS, Z, LDZ,
  312.      $                   IERR )
  313.             IF( IERR.GT.0 ) THEN
  314. *
  315. *              If DLAHQR failed to compute all NS eigenvalues, use the
  316. *              unconverged diagonal elements as the remaining shifts.
  317. *
  318.                DO 90 II = 1, IERR
  319.                   WR( I-NS+II ) = S( II, II )
  320.                   WI( I-NS+II ) = ZERO
  321.    90          CONTINUE
  322.             END IF
  323.          END IF
  324. *
  325. *        Form the first column of (G-w(1)) (G-w(2)) . . . (G-w(ns))
  326. *        where G is the Hessenberg submatrix H(L:I,L:I) and w is
  327. *        the vector of shifts (stored in WR and WI). The result is
  328. *        stored in the local array V.
  329. *
  330.          V( 1 ) = ONE
  331.          DO 100 II = 2, NS + 1
  332.             V( II ) = ZERO
  333.   100    CONTINUE
  334.          NV = 1
  335.          DO 120 J = I - NS + 1, I
  336.             IF( WI( J ).GE.ZERO ) THEN
  337.                IF( WI( J ).EQ.ZERO ) THEN
  338. *
  339. *                 real shift
  340. *
  341.                   CALL DCOPY( NV+1, V, 1, VV, 1 )
  342.                   CALL DGEMV( 'No transpose', NV+1, NV, ONE, H( L, L ),
  343.      $                        LDH, VV, 1, -WR( J ), V, 1 )
  344.                   NV = NV + 1
  345.                ELSE IF( WI( J ).GT.ZERO ) THEN
  346. *
  347. *                 complex conjugate pair of shifts
  348. *
  349.                   CALL DCOPY( NV+1, V, 1, VV, 1 )
  350.                   CALL DGEMV( 'No transpose', NV+1, NV, ONE, H( L, L ),
  351.      $                        LDH, V, 1, -TWO*WR( J ), VV, 1 )
  352.                   ITEMP = IDAMAX( NV+1, VV, 1 )
  353.                   TEMP = ONE / MAX( ABS( VV( ITEMP ) ), SMLNUM )
  354.                   CALL DSCAL( NV+1, TEMP, VV, 1 )
  355.                   ABSW = DLAPY2( WR( J ), WI( J ) )
  356.                   TEMP = ( TEMP*ABSW )*ABSW
  357.                   CALL DGEMV( 'No transpose', NV+2, NV+1, ONE,
  358.      $                        H( L, L ), LDH, VV, 1, TEMP, V, 1 )
  359.                   NV = NV + 2
  360.                END IF
  361. *
  362. *              Scale V(1:NV) so that max(abs(V(i))) = 1. If V is zero,
  363. *              reset it to the unit vector.
  364. *
  365.                ITEMP = IDAMAX( NV, V, 1 )
  366.                TEMP = ABS( V( ITEMP ) )
  367.                IF( TEMP.EQ.ZERO ) THEN
  368.                   V( 1 ) = ONE
  369.                   DO 110 II = 2, NV
  370.                      V( II ) = ZERO
  371.   110             CONTINUE
  372.                ELSE
  373.                   TEMP = MAX( TEMP, SMLNUM )
  374.                   CALL DSCAL( NV, ONE / TEMP, V, 1 )
  375.                END IF
  376.             END IF
  377.   120    CONTINUE
  378. *
  379. *        Multiple-shift QR step
  380. *
  381.          DO 140 K = L, I - 1
  382. *
  383. *           The first iteration of this loop determines a reflection G
  384. *           from the vector V and applies it from left and right to H,
  385. *           thus creating a nonzero bulge below the subdiagonal.
  386. *
  387. *           Each subsequent iteration determines a reflection G to
  388. *           restore the Hessenberg form in the (K-1)th column, and thus
  389. *           chases the bulge one step toward the bottom of the active
  390. *           submatrix. NR is the order of G.
  391. *
  392.             NR = MIN( NS+1, I-K+1 )
  393.             IF( K.GT.L )
  394.      $         CALL DCOPY( NR, H( K, K-1 ), 1, V, 1 )
  395.             CALL DLARFG( NR, V( 1 ), V( 2 ), 1, TAU )
  396.             IF( K.GT.L ) THEN
  397.                H( K, K-1 ) = V( 1 )
  398.                DO 130 II = K + 1, I
  399.                   H( II, K-1 ) = ZERO
  400.   130          CONTINUE
  401.             END IF
  402.             V( 1 ) = ONE
  403. *
  404. *           Apply G from the left to transform the rows of the matrix in
  405. *           columns K to I2.
  406. *
  407.             CALL DLARFX( 'Left', NR, I2-K+1, V, TAU, H( K, K ), LDH,
  408.      $                   WORK )
  409. *
  410. *           Apply G from the right to transform the columns of the
  411. *           matrix in rows I1 to min(K+NR,I).
  412. *
  413.             CALL DLARFX( 'Right', MIN( K+NR, I )-I1+1, NR, V, TAU,
  414.      $                   H( I1, K ), LDH, WORK )
  415. *
  416.             IF( WANTZ ) THEN
  417. *
  418. *              Accumulate transformations in the matrix Z
  419. *
  420.                CALL DLARFX( 'Right', NH, NR, V, TAU, Z( ILO, K ), LDZ,
  421.      $                      WORK )
  422.             END IF
  423.   140    CONTINUE
  424. *
  425.   150 CONTINUE
  426. *
  427. *     Failure to converge in remaining number of iterations
  428. *
  429.       INFO = I
  430.       RETURN
  431. *
  432.   160 CONTINUE
  433. *
  434. *     A submatrix of order <= MAXB in rows and columns L to I has split
  435. *     off. Use the double-shift QR algorithm to handle it.
  436. *
  437.       CALL DLAHQR( WANTT, WANTZ, N, L, I, H, LDH, WR, WI, ILO, IHI, Z,
  438.      $             LDZ, INFO )
  439.       IF( INFO.GT.0 )
  440.      $   RETURN
  441. *
  442. *     Decrement number of remaining iterations, and return to start of
  443. *     the main loop with a new value of I.
  444. *
  445.       ITN = ITN - ITS
  446.       I = L - 1
  447.       GO TO 50
  448. *
  449.   170 CONTINUE
  450.       RETURN
  451. *
  452. *     End of DHSEQR
  453. *
  454.       END
  455.